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ABSTRACT 


The  Arctic  Ocean  has  been  a  subject  of  increasing  interest 
in  recent  years  due  to  the  reduction  of  the  sea-ice 
thickness  and  spatial  coverage  and  its  implications  for 
climate  change.  The  future  state  of  the  Arctic  is  likely 
to  be  linked  to  vertical  heat  transport  by  microscale 
processes,  specifically,  double-diffusive  convection.  A 
series  of  realistic  three-dimensional  direct  numerical 
simulations  (DNS)  were  conducted  to  assess  the  vertical 
heat  transport  through  thermohaline  staircases  in  the 
Arctic  region.  Results  revealed  that  vertical  fluxes 
exceeded  those  of  extant  "four-thirds  flux  laws"  by  as  much 
as  a  factor  of  two,  and  suggest  that  the  4/3  exponent 
requires  downward  revision.  Results  also  showed  that  two- 
dimensional  DNS  can  provide  an  accurate  approximation  of 
heat  fluxes  when  the  density  ratio  is  sufficiently  large. 
DNS  results  also  reveal  that  the  models  with  rigid 
boundaries  result  in  heat  flux  estimates  that  are  lower 
than  those  from  models  with  periodic  boundary  conditions. 
Finally,  the  DNS-derived  flux  law  was  applied  to  Arctic 
data  and  results  supported  the  conclusion  that  lab-derived 
flux  laws  significantly  underestimate  heat  flux.  All  of 
these  results  suggest  that  vertical  heat  transport  due  to 
double-diffusive  convection  is  a  significant  contributor  to 
the  measured  reduction  of  Arctic  sea-ice. 
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I.  INTRODUCTION  AND  BACKGROUND 

A.  ARCTIC  OCEAN  AND  CLIMATE  CHANGE 

Recent  dramatic  changes  in  Arctic  sea-ice 
characteristics,  specifically  its  extent  and  thickness, 
have  sparked  a  wide  range  of  scientific  interest  and 
concern.  There  has  been  greater  focus  on  this  fiercely 
debated  subject  as  research  activities  increase  in  response 
to  the  reduction  in  sea-ice  coverage.  These  activities 
reflect  not  only  the  scientific,  but  also  the  socioeconomic 
importance  of  a  reduced  Arctic  sea-ice  region,  with 
increased  interest  in  a  wide  range  of  disciplines  including 
merchant  shipping,  maritime  safety  and  security,  and 
resource  management.  Figure  1  illustrates  the  dramatic 
declination  of  sea-ice  extent,  especially  in  the  late 
summer  of  2012  (the  dashed  line) .  Combined  estimates  of  ice 
thickness  derived  from  submarine  records,  as  well  as  ice, 
cloud  and  land  elevation  satellite  laser  altimetry,  provide 
a  long-term  record  of  sea-ice  thickness;  these  data  suggest 
a  decrease  in  the  mean  overall  thickness  of  the  sea-ice 
over  a  region  covering  about  38%  of  the  Arctic  Ocean 
(Gardner  et  al .  2012)  . 

Observations  suggest  that  there  is  a  thermodynamic 
coupling  between  heat  in  the  ocean  interior  and  sea-ice 
(Polyakov  et  al .  2010)  .  This  coupling  is  supported  by 
arguments  for  the  role  of  oceanic  heat  in  shaping  the 
Arctic  sea-ice,  linked  to  the  presence  of  warm 
intermediate-depth  water  of  Atlantic  origin.  "Atlantic 
Water"  (AW)  carries  a  vast  amount  of  heat  that  observations 
have  shown  dissipates  away  from  its  Arctic  entrance,  along 
various  spreading  pathways. 
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Arctic  Sea  Ice  Extent 

(Area  of  Ocean  with  at  least  15)6  sta  ice) 


Figure  1.  Arctic  sea-ice  extent  declination.  Note  2012 
data  make  up  the  dashed  green  line  (From 
Fetterer  et  al .  2012)  . 


While  some  AW  heat  is  lost  along  spreading  pathways, 

some  of  this  heat  must  be  lost  to  the  overlying  halocline 

waters  (Rudels  et  al .  1996;  Steele  and  Boyd  1998;  Martinson 

and  Steele  2001)  .  Previous  analytical  estimates  of  upward 

heat  flux  from  the  AW  yield  values  ranging  from  1. 3-2.1 

W/m2,  while  parameterization  for  the  bulk  heat  transfer 

coefficient  from  the  mixed  layer  temperature  has  yielded 

vertical  heat  flux  of  3-4  W/m2  (Krishfield  and  Perovich 

2005)  .  The  diversity  of  these  oceanic  heat  flux  estimates 

can  be  attributed  to  strong  high-latitude  temporal 
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variability,  the  various  warm  and  cold  periods  experienced 
throughout  the  twentieth  century,  and  the  continually 
developing  understanding  of  various  heat  advection 
mechanisms.  AW  warms  and  cools  in  response  to  different 
climate  periods,  and  the  remnants  of  these  signals  arrive 
and  are  felt  in  the  Arctic  Ocean  years  later.  Observations 
from  the  2000s  document  a  large-scale  warming  of  the  AW 
layer  with  temperatures  unprecedented  in  the  history  of 
regional  instrumental  observations  (Polyakov  et  al .  2010). 


0.3  1.0  2.0  -1  -0.5  0  0.5  1 

Figure  2.  (a)  Mean  AW  temperature  (°C)  averaged  over 

the  1970s;  (b)  AW  temperature  anomalies  (°C) 

averaged  over  the  1990s,  computed  relative 
to  the  climatological  mean;  (c)  same  as  (b) 
but  for  2007  data  only  (From  Polyakov  et  al . 
2010) . 


Using  climatology.  Figure  2  clearly  depicts  the 
anomalously  warm  AW  for  the  1990s  and  more  recently,  for 
2007.  At  the  time,  2007  held  the  record  for  minimum  summer 
sea-ice  extent  in  the  Arctic  region.  An  anomalous  warming 
event  similar  to  the  one  documented  in  2007  would  likely  be 
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observed  in  2012;  in  September  2012,  Arctic  sea-ice  reached 
a  new  minimum  extent  (Figure  1) . 

It  was  observed  that  sea-ice  thickness  was  also 
shrinking.  Regions  of  the  Beaufort  and  Laptev  Seas  saw  ice- 
free  areas  opening  faster  than  before,  indicative  of  the 
thinner  ice  melting.  Observations  from  NASA's  Operation 
IceBridge  (collection  of  airborne  remote  sensing 
measurements)  were  used  to  derive  sea-ice  thickness  and 
revealed  thinning  in  the  Beaufort  Sea  that  could  be 
attributed  to  warmer  AW  (Figure  3) . 


0  1  2  3  4  S 

Sea  Ice  Thickness 


Figure  3.  Arctic  sea-ice  thickness  derived  from  2012 
Operation  IceBridge  (NSIDC  website  2012). 

Remote  sensing  observations  also  reveal  that  the 
thickest  (oldest)  Arctic  sea-ice  is  thinning  and  melting  at 
a  faster  rate  than  the  younger  (thinner)  ice  at  the  edges 
of  the  Arctic  Ocean's  floating  ice  cap.  Multi-year  ice 
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typically  survives  the  cyclical  summer  melt  season,  but 
this  older  ice  has  diminished  with  each  passing  winter  over 
the  past  three  decades  (Figure  4) .  The  average  thickness  of 
the  Arctic  sea-ice  cover  is  declining  because  it  is  rapidly 
losing  its  thick  component.  These  dramatic  changes  point  to 
the  possibility  of  a  warmer  ocean  interior  impacting  the 
melting  rate  of  sea-ice. 


Figure  4 .  NASA  remote  sensors  show  the  thickest  parts 
of  the  Arctic  Ice  Cap  are  shrinking  (After 
NASA  website  2012)  . 


1 .  Heat  Required  to  Melt  Arctic  Sea-Ice 

The  observed  trend  of  a  warming  AW  layer  in  the  Arctic 
region,  coupled  with  the  irrefutable  evidence  of  shrinking 
Arctic  sea-ice  coverage,  elicits  the  question:  How 
significant  do  these  AW  temperature  changes  have  to  be  to 
impact  the  upward  flux  of  heat  and,  in  turn,  the  sea-ice? 

Additionally,  how  much  heat  is  needed  to  melt  the  sea-ice 

in  the  Arctic  Ocean,  and  could  it  be  supplied  by  the 
intruding  AW?  The  heat  required  to  melt  ice  one  meter  thick 

is  8  kcal/cm2  (3.35xl08  J/m2)  for  a  latent  heat  of  fusion  of 

80  cal/gm  (Turner  2010)  .  Aagaard  and  Greisman  (1975)  have 
calculated  the  average  net  rate  of  transport  of  heat  into 
the  Arctic  by  the  Atlantic  inflow  to  be  about  6.7xl010  kW, 
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with  a  mean  temperature  of  2.2  °C .  This  current  accounted 
for  more  than  half  the  total  heat  flux  through  all  passages 
into  the  Arctic  during  a  yearlong  period  (1971-1972) .  More 
recent  observations  find  that  Aagaard  and  Greisman' s  result 
underestimates  subsurface  advective  heat  flux  since  it  does 
not  factor  in  the  more  recent  increase  in  temperature  of 
the  AW  and  its  extension  through  the  Arctic  (Turner  2010) . 

If  the  above-referenced  heat  flux  of  6.7xl010  kW  is 
distributed  uniformly  across  the  entire  Arctic  Ocean,  which 
has  a  surface  area  of  approximately  107  km2,  and  the  average 
heat  loss  from  this  layer  is  about  6.7  W/m2,  then  it  can  be 
concluded  that  in  12-18  months  there  is  enough  heat  in  the 
AW  alone  to  melt  one  meter  of  sea  ice.  Furthermore,  there 
is  enough  heat  from  the  intruding  AW  to  melt  all  the  sea 
ice  with  a  mean  thickness  of  2.5  m  in  a  period  of  four 
years  (Turner  2010)  .  This  assumption  is  only  feasible  if 
there  is  a  mechanism  for  transferring  the  heat  from  the 
intermediate  water  level  to  the  surface  under  the  ice  pack. 
These  mechanisms  that  could  bring  heat  from  the  subsurface 
layer  to  the  surface  over  this  very  short  time  scale  are 
the  subject  of  several  studies  (Turner  and  Veronis  2004) . 

2 .  Mechanisms  that  Impact  Melting  Rate  of  Sea  Ice 

Kwok  and  Untersteiner  (2011)  noted  that  the  surplus 
heat  needed  to  explain  the  complete  loss  of  Arctic  sea-ice 
during  the  past  few  decades  is  on  the  order  of  1  W/m".  As 
the  Arctic  sea-ice  continues  to  melt  at  an  accelerated 
rate,  the  mechanisms  responsible  for  this  phenomenon  are  of 
increasing  interest.  Perhaps  the  most  extensively  studied 
mechanism,  especially  as  a  target  for  determining  evidence 

of  climate  change,  is  the  heat  input  from  the  atmosphere 
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directly  to  the  surface  of  the  ice  and  the  surface  mixed 
layer.  Observations  gathered  during  the  year-long  Surface 
Heat  Budget  of  the  Arctic  Ocean  (SHEBA)  study  revealed  an 
annual  average  net  heat  flux  into  undeformed  ice  of  7.5  W/m2 
(Perovich  et  al .  2003)  .  The  SHEBA  values  are  2-3  times 

larger  than  prior  results,  which  were  typically  around  3 
W/m2.  This  increase  in  atmospheric  input  can  help  explain 
the  increased  melting  rate  observed  in  recent  years. 

The  fact  that  there  is  enough  heat  in  the  AW  to  melt 
all  the  sea-ice  in  the  Arctic  Ocean  in  just  a  few  years 
points  to  the  importance  of  oceanic  advective  processes 
that  can  bring  this  heat  more  rapidly  to  the  surface  mixed 
layer.  One  such  mechanism  is  the  mixing  caused  by  the  down¬ 
flow  of  dense  shelf  water.  A  two  dimensional  (2D)  density 
current  flowing  into  a  two-layer  stratified  basin  will 
follow  along  the  downward  slope,  causing  a  dilution  due  to 
the  entrainment  of  lighter  fluid  above.  In  the  Arctic 
context,  the  rejection  of  brine  on  the  shelves  during 
winter  could  lead  to  a  dense  layer  spreading  out  along  the 
halocline,  increasing  in  density  and  eventually  penetrating 
below  the  level  of  the  core  of  the  Atlantic  water.  This 
could  eventually  bring  the  warmer  AW  into  contact  with  the 
surface  sea-ice  through  a  combination  of  upwelling  and  the 
erosion  of  the  existing  stratification  that  allows  the 
surface  forcing  to  penetrate  more  deeply  into  the  water 
column  (Turner  2010)  . 

Another  mechanism  that  impacts  the  rate  of  ice-melt 
was  the  subject  of  extensive  experimentation  by  Turner  and 
Veronis  (2004)  .  They  examined  the  effect  of  the  vertical 
transport  of  heat  from  warmer  layers,  intruding  laterally 
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below  the  surface  into  a  salinity  gradient,  on  the  melting 
sea-ice  floating  on  the  surface.  The  melting  rate  was  found 
to  be  dependent  on  the  magnitude  of  the  salinity  gradient 
and  on  the  presence  (or  absence)  of  heating  from  below.  The 
melting  rate  decreased  (increased)  as  the  salinity  gradient 
was  increased  (decreased)  .  These  results  demonstrate  the 
potential  effect  of  the  vertical  double  diffusive  transport 
of  heat,  the  mechanism  of  focus  in  this  thesis. 

While  solar  radiation  is  undoubtedly  the  biggest 
contributor  to  the  reduction  of  Arctic  sea-ice, 
oceanographic  advection  processes  likely  play  a  much  larger 
role  than  previously  believed,  thus  explaining  the  current 
heat  surplus  in  the  Arctic  and  illustrating  the  larger  role 
diffusive  convection  plays  in  transporting  AW  vertically  to 
the  surface.  Figure  5  depicts  the  interaction  of  the  three 
primary  mechanisms  discussed  in  this  section,  and  the 
current  estimates  of  heat  flux,  based  on  previous  studies. 
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Figure  5.  1:  Upward  double  diffusion  transport  from 

the  AW.  2:  Heat  flux  from  the  atmosphere 
into  the  surface  mixed  layer.  3:  Mean 
surface  heat  flux  over  the  whole  region 
required  to  balance  the  net  input  of  heat  by 
the  intruding  AW  (After  Turner  2010)  . 


B.  DOUBLE -DIFFUSIVE  CONVECTION 

Turner  (1973)  describes  double-diffusive  convection  as 
a  phenomenon  that  occurs  in  a  fluid  with  two  (or  more) 
components  with  different  molecular  dif fusivities .  With 
regard  to  the  ocean,  net  differences  in  the  diffusivity  of 
heat  and  salt  are  of  interest,  and  give  rise  to  two 
regimes:  salt  fingering  and  diffusive  convection. 

1 .  Salt  Finger  Regime 

Salt  fingering  instability  is  prevalent  in  low  and 
mid-latitudes,  where  the  combination  of  evaporation  and 
radiative  heating  at  the  ocean  surface  results  in  a  layer 
of  warm,  salty  water  overlying  a  comparatively  cooler. 
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fresher  layer  of  water.  These  layers  are  thermally  stably 
stratified,  but  have  an  inhomogeneous  (destabilizing) 
composition.  Since  heat  diffuses  one  hundred  times  faster 
than  salt,  fluid  displaced  downward  rapidly  loses  its  heat 
excess  while  maintaining  its  larger  salt  concentration.  It 
becomes  denser  than  the  surrounding  environment  and  continues 
to  sink,  forming  salt-finger  structures  (Figure  6) . 


Figure  6.  Salt  fingers.  Warm,  salty  layer  atop  a  cool, 
fresh  layer  results  in  fingers  of  salty, 
cooler  water  extending  downwards  from  the 
layer  interface  (Image  from  Garaud  Research 
website  2012 )  . 

2 .  Diffusive  Convection  Regime 

In  the  second  case,  diffusive  instability,  the  fluid 
is  stably-stratified  despite  an  unstable  thermal  gradient. 
It  is  prevalent  in  high-latitude  locations  where  cold  and 
fresh  water  overlies  relatively  warm  and  salty  water.  This 
phenomenon  can  occur  in  the  ocean  under  the  ice-shelf, 
where  the  slow  melting  of  the  ice  releases  cool  fresh  water 
over  the  warmer,  saltier  water  below.  This  thesis  is 
focused  on  the  diffusive  convection,  commonly  observed  in 
the  Beaufort  Gyre  region  of  the  Arctic  Ocean. 
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To  understand  the  dynamics  of  diffusive  convection, 
consider  a  parcel  of  warm,  salty  water  that  is  displaced 
upward  across  a  thin  interface  into  cooler,  fresher  water. 
The  parcel  quickly  diffuses  its  heat  into  the  new 
surroundings  but  slowly  releases  its  salt  content.  The 
parcel  becomes  denser,  so  it  sinks  back  across  the 
interface  and  overshoots  its  original  starting  point.  Now 
the  parcel  has  returned  to  warmer  waters,  and  will  start 
this  entire  process  all  over  again.  This  over-stable 
process,  known  as  oscillatory  diffusive  convection, 
transports  heat  upward  in  an  oscillatory  fashion  (See 
Figure  7 ) . 

This  mixing  results  in  the  formation  of  layers  of  near 
constant  temperature  and  salinity— diffusive  layering.  Well- 
mixed  layers  separated  by  sharp  diffusive  interfaces  are 
responsible  for  the  stair-like  appearance  of  temperature 
profiles  from  the  central  Arctic  thermocline.  Such  profiles 
are  referred  to  as  thermohaline  staircases  and  this  process 
is  known  to  be  prevalent  in  the  Beaufort  Sea  region  of  the 
Arctic  Ocean.  The  physical  processes  associated  with  the 
maintenance  of  this  structure  are  as  follows:  Heat  and  salt 
are  transferred  upward  through  the  diffusive  interface.  The 
water  above  the  interface  gains  enough  heat  to  begin  upward 
motion  due  to  excess  buoyancy.  This  water  begins  to  rise 
within  the  layer,  which  drives  the  convection  that 
ultimately  maintains  the  homogeneous  properties  of  the 
mixed  layer.  At  the  same  time,  the  particles  just  below  the 
interface  cool  off  and  sink,  which  also  drives  the 
convection  cell  in  that  layer.  This  process  maintains  a 
homogeneous  distribution  of  T-S  in  each  layer  (Figure  8)  . 
The  vertical  transport  of  heat  via  diffusive  layering  is 
the  key  mechanism  of  interest  being  studied  in  this  thesis. 
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C.  FLUX  LAWS  AND  HEAT  FLUX  ESTIMATES 

For  systems  characterized  by  active  diffusive 
layering,  the  interaction  between  layers  is  controlled  by 
the  salt  and  heat  fluxes  across  the  diffusive  interfaces— a 
dependency  that  has  motivated  extensive  study.  After 
conducting  laboratory  experiments.  Turner  (1965)  proposed 
that  the  balance  of  temperature  variation  between  layers 
and  heat  flux  is  given  by 

4 

aFT  =  A,  (orAr ( 1 ) 


„  J3AS 

RP=~ 77  (2) 

aAT 

Based  on  the  exponent.  Equation  (1)  is  commonly 
referred  to  as  the  "four-thirds  flux  law,  "  where  FT  is 
temperature  flux,  a  is  the  thermal  expansion  coefficient,  (3 
is  the  saline  contraction  coefficient,  and  Ai  has  the 
dimensions  of  velocity  and  is  a  function  of  the  density 
ratio  (Rp)  .  AT  and  AS  are  the  differences  in  temperature  and 
salinity  between  adjacent  layers.  This  law  states  that  the 
fluxes  associated  with  double  diffusion  are  controlled 
entirely  by  the  variation  in  temperature  or  salinity  across 
the  diffusive  interface  and  is  independent  of  layer 
thickness.  Using  dimensional  analysis,  the  4/3  flux  law  can 
be  written  as 


ccFT=C{Rp) 


■(ccATy 


where  the  (unknown)  coefficient  C  (Rp)  is  a  function  of 
density  ratio  only,  g=9.8  m/s,  v=1.8xl0~6  m2/s  is  the 
kinematic  viscosity,  and  k=1.4x1CT7  m2/s  is  the  molecular 
diffusivity  of  heat. 
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In  1976,  Marmorino  and  Caldwell  examined  Turner's 
four-thirds  flux  law  in  more  depth  by  investigating  a 
larger  range  of  Rp  values  and  additional  initial  conditions. 
They  also  factored  in  Huppert' s  (1971)  heat  flux  law  on 
interface  stability 

H/Hsp=3.8Rp~2  (4) 

and  heat  flow 

i  1 

Hsp  =  (0.085) kT(g a  /  AT3  (5) 

through  a  fictitious  non-def ormable  interface  (Turner 

1965)  .  Using  Huppert' s  formulae,  they  determined  the 
relationship  between  H  and  Hsp  to  be 

—  =  0.101exp{4.6exp[-0.54(/^  -l)]J  (6) 

Hsp 

and  proposed  the  new  empirical  formulation  for  the 
coefficient  C  (Rp) 

C(/?p)  =  0.00859exp(4.6exp[-0.54{/?p-l}])  (7) 

Using  empirical  formulations,  Kelley  (1990)  found  that 
at  low  density  ratios  (Rp<2.5),  temperature  fluxes  were 
significantly  smaller  than  those  given  by  Marmorino  and 
Caldwell  (1976)  .  This  lead  to  Kelley  proposing  the 
formulation 

C(^)  =  0.0032exp(4.8/«y)  (g) 

based  on  what  he  referred  to  as  a  full  collection  of 
laboratory  experiments;  he  analyzed  experiments  for  1<RP<10, 
while  Turner  (1973)  had  limited  his  study  to  Rpd7 . 


Kelley 

also  suggested 

that 

the 

exponent 

to 

use 

in 

Equation  (1) 

i  be  reduced 

to 

5/4. 

However,  (3) 

with 

either 

(7)  or  (8) 

continue  to 

be 

used 
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as  a 

standard 

method 

for 

calculating  fluxes.  Additional  analysis  is  required  to 
support  or  challenge  the  existing  four-thirds  flux  law. 

D .  RESEARCH  MOTIVATION 

This  study  is  focused  on  the  diffusive  convection 
regime  of  double  diffusion.  The  area  of  research  utilized 
high  resolution  three  dimensional  (3D)  numerical  modeling 
aimed  at  determining  the  vertical  heat  transport  through 
double  diffusive  staircases  in  the  western  Arctic  Ocean. 
Direct  numerical  simulations  (DNS)  guided  and  validated 
attempts  to  formulate  more  realistic  and  accurate  heat 
flux,  anticipating  that  laboratory  derived  flux  laws  do  not 
adequately  represent  the  vertical  heat  transport  in  the 
Arctic  Ocean.  This  research  determined  whether  it  is 
possible  to  reliably  assess  the  performance  of  flux  models 
using  DNS,  especially  for  experiment  sizes  that  most 
accurately  reflect  an  actual  thermohaline  step  height. 

The  model  runs  were  high  resolution  3D  experiments, 
which  were  computationally  demanding.  This  higher  level  of 
resolution  helped  determine  how  accurate  earlier  2D 
numerical  models  are  in  terms  of  predicting  vertical  heat 
flux.  If,  for  example,  there  was  not  a  significant 
difference  in  the  results  between  high  resolution  3D 
modeling  and  2D  modeling,  fellow  researchers  would  not 
worry  about  acquiring  an  extensive  amount  of  computational 
resources.  On  the  other  hand,  if  the  results  were 
significantly  different,  then  clearly  2D  experimentation  is 
not  an  adequate  approach  to  modeling  this  microscale 
feature . 
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Another  element  examined  in  this  study  was  the  role  of 
boundary  conditions  in  controlling  the  heat  transport  in 
staircase  models.  For  example,  what  form  of  boundary 
conditions  best  represents  the  dynamics  of  Arctic 
staircases?  This  was  simulated  by  adding  boundaries  to  the 
DNS  that  were  designed  to  replicate  the  restrictions 
encountered  in  laboratory  experiments.  Comparing  the 
results  of  DNS  with  periodic  and  rigid  boundaries 
demonstrated  the  sensitivity  of  modeled  fluxes  to  boundary 
conditions.  Finally,  the  modeled  results  were  validated  by 
incorporating  real  input  values  taken  from  Ice-Tethered 
Profiler  (ITP)  data  to  determine  applicability  of  DNS- 
derived  flux  laws  in  the  ocean. 

Of  note,  this  is  the  first  series  of  fully  3D  DNS  to 
be  performed  for  the  diffusive  convection  for  the  step-size 
and  diffusivity  ratio  that  is  representative  of  Arctic 
staircases.  It  was  anticipated  that  DNS-derived  flux  values 
would  be  significantly  higher  than  those  previously 
encountered  in  laboratory  experiments,  and  that  it  is 
feasible  for  vertical  heat  transport  via  double  diffusion 
to  be  significant  enough  to  impact  the  melting  of  Arctic 
sea-ice . 
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II.  MODEL  DESCRIPTION  AND  SETUP 


DNS  were  conducted  in  both  two  and  three  dimensions. 
Lower  resolution,  2D  experiments  were  conducted  on  the 
Naval  Postgraduate  School  (NPS)  high-performance  Hamming 
Linux  cluster  using  between  64  and  256  processors  per  run. 
3D  model  runs  were  performed  on  the  Department  of  Defense's 
High  Performance  Computing  Modernization  Program  (HPCMP) 
Open  Research  Systems,  Cray  XE6  (Chugach)  and  Cray  XE6m 
(Copper)  located  at  the  Engineer  Research  and  Development 
Center  (ERDC)  in  Vicksburg,  MS;  and  on  the  University  of 
Texas  at  Austin' s  Advance  Computing  Center  (TACC)  using 
their  Sun  Constellation  Linux  Cluster  (Ranger) ,  and  their 
Dell  Linux  Cluster  (Lonestar) .  The  resources  used  varied 
from  256-2048  processors  per  run  for  computing  and  model 
efficiency  testing;  results  determined  that  512  processors 
are  optimal  for  final  high  resolution  runs. 

The  simulations  were  performed  using  the  Parallel 
Algorithm  for  Double  Diffusive  Instability  (PADDI)  model,  a 
fully  de-aliased,  pseudo-spectral,  Fortran  numerical  code 
developed  specifically  for  double  diffusive  convection  and 
which  utilizes  a  domain  that  is  periodic  in  each  spatial 
dimension  (Lx,  Ly,  Lz)  .  The  model  was  parallelized  using  the 
Message  Passing  Interface  (MPI),  and  code  was  compiled 
using  the  Intel  Fortran  compiler.  Experiments  were  run  in 
single  precision  once  it  was  verified  that  the  differences 
between  single  and  double  precision  results  were  not 
significant.  The  2D  study  examined  diffusive  convection  in 
the  horizontal  (x)  and  vertical  (z)  directions. 
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For  modeling  purposes,  temperature  and  salinity  output 
fields  were  decomposed  into  the  linear  basic  state  ( T,S )  and 
perturbations  (T,S) .  The  relevant  equations  are  then  the 
Boussinesq  equations  of  motion,  expressed  in  terms  of  T  and 
S  and  non-dimensionalized  using  the  "standard"  double- 
diffusive  system  (Radko  2013)  : 

/  =  (J^V)4  (9) 

gaZ 

kt 

V  =  y-  (10) 

l2 

t  =  ~  (ID 

KT 

Here,  1  denotes  the  length  scale,  v  is  the  velocity  scale, 
t  is  the  time  scale,  Tz  is  the  temperature  gradient  and  aTzl 
is  the  scale  for  both  temperature  and  salinity 
perturbations  (Radko  and  Stern  1999).  These  relationships 
lead  to  the  following  governing  non-dimensional  equations 
that  are  solved  by  PADDI 


✓^11 

—  +  u-Vu)  =  -Vp  +  (T-lS)k  +  V2u 
dt 

(12 

V-u  =  0 

(13 

— -+u-vr-ku  =  v2r 

dt 

(14 

dS 

— +  u-V5-i?1ku=rV25 

dt  p 

(15 

for  initial  conditions  representing  the  system  at  rest 

u(r,  t  =  0)  (16) 
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where  the  step-like  distribution  of  temperature  and 
salinity  are  the  result  of  slight  perturbations  caused  by 
random  computer-generated  noise. 


Equation  (12)  is  the  equation  of  state.  Equation  (13) 
is  the  continuity  equation,  and  Equations  (14)  and  (15)  are 
the  thermal  and  compositional  advection-dif fusion 
equations.  Here,  u  denotes  the  velocity  field  and  p  is 
pressure  (Peyret  2002;  Canuto  et  al .  2007;  Stellmach  and 

Hansen  2008) .  Equations  (12) -(15)  can  be  solved  in  either  a 
2D  or  3D  periodic  domain,  thereby  representing  conditions 
in  the  real  ocean  where  the  influence  of  boundary 
conditions  are  minimal.  The  governing  parameters  are  the 
Prandtl  number  (Pr) ,  molecular  diffusivity  ratio  (t),  and 
Rp.  Pr  and  i  are  defined  as 


kt 


(17) 


where  kt  is  the  molecular 
molecular  diffusivity  of 


dif fusivi 
salt . 


ty 


of 


heat 


(18) 

and  ks  is  the 


One-layer  DNS  were  conducted  using  parameter  settings 
that  reflect  those  found  in  the  vicinity  of  the  Beaufort 
Gyre  in  the  Arctic  Ocean,  where  recently,  ITPs  have 
recorded  temperature  and  salinity  profiles  that  clearly 
show  stepped  patterns  at  depths  200-300  m  (Timmermans  et 
al .  2008).  Therefore,  for  all  experiments,  v=1.8xl0~6  m2/s, 

kt=1.8><10~6  m2/s,  ks=7x10~10  m2/s,  yielding  i=0.005  and  Pr=13. 

Previous  DNS  were  unable  to  use  such  a  small  i  value 
because  of  the  computational  difficulty  of  resolving  the 
extremely  fine  salt  dissipation  scale  (e.g.  Prikasky  2007; 
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Caro  2009;  Carpenter  et  al .  2012)  .  From  ITP  data,  the 
average  values  a=5.625xl0~s  1/°C  and  T_  =1 . 3  61  x  1 0~2  °C/m  were 

determined  over  the  appropriate  depth  range.  Thus, 
l=1.35xl0~2  m  and  the  relevant  non-dimensional  layer¬ 
thickness  range  is  74-740.  The  density  ratio  in  the 
diffusive  layer  of  the  Arctic  Ocean  is  typically  in  the 
range  2hRp<lO  and  this  analysis  is  based  on  a  series  of 
experiments  with  Rp= [2, 3, 4, 6, 8, 10]  . 

For  2D  experiments,  the  chosen  non-dimensional  domain 
size  was  300x300  so  as  to  reflect  step-sizes  in  an  actual 
Arctic  staircase.  The  corresponding  mesh  size  was  756x3072 
to  sufficiently  resolve  the  scale  of  salt  dissipation  for 
the  realistic  diffusivity  ratio  used.  For  3D  experiments, 
the  typical  non-dimensional  domain  size  was  300x150X300  and 
the  corresponding  mesh  size  was  756x384x3072.  The  box  size, 
H,  is  varied  in  one  set  of  2D  and  3D  experiments. 

One  further  experiment,  designed  to  replicate  imposed 
laboratory-type  boundary  conditions  and  observe  their 
influence  on  vertical  heat  transport  in  thermohaline 
staircases,  was  conducted.  For  this  experiment,  rigid 
stress-free  boundaries  were  assumed  in  all  spatial 
dimensions,  with  minimal  changes  to  the  underlying  code. 
Table  1  outlines  all  the  DNS  conducted  in  this  study. 
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2/3D 

Type 

Size  (m) 

Mesh  Size 

T 

Rp 

Total  #  of 
experiments 

3D 

Periodic 

4x2x4 

768x384  x3072 

.005 

2,3,4,6.8.10 

6 

3D 

Periodic 

4x2x2 

768x384x1536 

.005 

4 

1 

3D 

Periodic 

4x2x1 

768x384x768 

.005 

4 

1 

3D 

Periodic 

4x2x0. 5 

768x384x384 

.005 

4 

1 

2D 

Periodic 

4x4 

768x3072 

.005 

2,3,4,6,8,10 

6 

2D 

Periodic 

4x2 

768x1536 

.005 

4 

1 

2D 

Periodic 

4x1 

768x768 

.005 

4 

1 

2D 

Periodic 

4x0.5 

768x384 

.005 

4 

1 

2D 

Non¬ 

periodic 

4x4 

768x3072 

.005 

1.5+ 

1 

Table  1.  A  list  of  DNS  conducted  in  this  study. 


Prior  to  this  study,  the  computational  demands  of 
performing  DNS  in  3D  led  modelers  to  either  resort  to  2D 
simulations  or  to  compromise  in  their  choice  of  parameters 
(use  unrealistic  values  of  diffusivity  ratio  or  step 
heights,  for  instance)  .  Caro  (2009)  performed  DNS  of  a 

fluid  with  periodic  boundary  conditions  but  with 

unrealistic  i  and  H  values  (relative  to  what  is  typically 
observed  in  the  Arctic  Ocean)  because  he  found  using 

smaller  diffusivity  ratios  and  larger  grid  sizes  to  be 
computationally  prohibitive.  Because  of  these  limitations, 
he  extrapolated  high-tau  simulations  to  evaluate  the 
vertical  heat  transport  for  Arctic  conditions  and  still 

found  it  to  be  significantly  larger  than  inferred  from 
laboratory  experiments.  He  concluded  that  laboratory 
experiments  being  conducted  in  small  tanks  that  have  rigid 
boundaries  would  have  guite  different  parameter  ranges  and 
were  very  susceptible  to  experimental  errors. 
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Carpenter  et  al .  (2012)  also  performed  a  series  of 

three-dimensional  DNS,  but  at  a  significantly  higher 
resolution,  utilizing  a  larger  domain  (Lx,  Ly,  Lz)  and  mesh 
size,  and  a  more  realistic  range  of  i  .  Their  results  also 
showed  that  the  flux  laws  of  Kelley  (1990)  gave  lower 
predicted  fluxes  than  those  observed  in  their  simulations, 
with  a  stronger  dependence  on  density  ratio.  Their 
experiments  also  had  imposed  unrealistic  rigid  boundary 
conditions  that  might  have  stunted  the  flux  values.  Both 
Caro's  (2009)  and  Carpenter  et  al . ' s  (2012)  results  bring 
into  question  the  relevance  of  lab-based  flux  laws  when 
applied  to  Arctic  staircases. 

To  that  end,  the  key  to  this  study  was  to  incorporate 
all  the  relevant  parameters  and  measures  that  would  most 
accurately  represent  a  typical  Arctic  thermohaline  step.  To 
date,  no  3D  DNS  has  simultaneously  utilized  a  mesh  size 
that  sufficiently  resolves  the  scale  of  salt  dissipation 
for  a  realistic  t=.005;  increased  domain  size  to  accurately 
reflect  step-sizes  in  an  actual  staircase  (H~4  m)  ;  allowed 
the  DNS  to  run  for  a  time  period  that  ensures  accurate 


diagnostics 

;  and 

guaranteed 

that 

all 

parameters 

and 

conditions , 

such 

as  a  lack 

of 

rigid 

boundaries , 

are 

reflective 

of  those  found  in 

the 

Arctic 

Such  a  set 

-up 

required  an  unprecedented  amount  of  resources  but  produced 
flux  law  derivations  that  are  more  representative  of  real- 
world  conditions. 
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Ill .  MODEL  RESULTS 


A.  TWO-DIMENSIONAL  MODEL 

1.  General  Characteristics 

Modeling  in  2D  served  three  purposes.  First,  it 
provided  a  testing  platform  for  checking  all  the 
parameters,  and  ensured  the  model  output  met  expectations 
prior  to  running  the  eguivalent  3D  jobs  (e.g.  flux  trends 
and  dependencies) .  As  anticipated,  computational  resources 
needed  for  running  high  resolution  3D  jobs  were  limited  and 
could  not  be  wasted.  Second,  2D  model  runs  were  conducted 
so  that  the  results  could  be  compared  to  their  3D 
counterparts  to  determine  the  dependencies  of  heat  flux  on 
step  heights  (H)  and  Rp.  Finally,  the  last,  rigid  boundary 
model  run  was  completed  in  2D  in  favor  of  its  numerical 
efficiency.  Refer  to  Table  1  for  the  complete  details  of 
the  10  experiments  ran  in  2D. 

Figures  9-11  show  a  full  temperature  field  plot 
revealing  the  spatial  structure  of  the  diffusive 
convection.  The  two  layered  fluid  is  initially  at  rest 
(Figure  9),  with  no  activity  across  the  single  interface. 
Higher  temperature  values  are  shown  in  red  and  lower  values 
in  blue.  Figure  10  depicts  signs  of  initial  convective 
activity  across  the  interface  as  the  heat  exchange  begins 
to  diffuse  between  the  two  layers.  Figure  11  more  clearly 
shows  the  heat  signatures  of  active  transport  and  the 
temperature  field  becoming  better  mixed  just  above  and 
below  the  interface.  In  Figure  11,  the  field  is  "zoomed  in" 
on  a  particular  area  of  the  interface  to  better  show  the 
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detail  of  enhanced  convective  dynamics,  further  confirmed 
by  plumes  emanating  from  the  interface. 


Figure  9.  2D  temperature  field  at  t=0 .  Simulation  is 
at  rest  and  no  convection  has  begun.  Plot 
aligns  with  the  x-z  axis  (not  to  scale) . 


Figure  10.  2D  temperature  perturbation  field  at  t=27 
minutes.  Convective  plumes  become  visible  in 
both  layers. 


Figure  11.  2D  temperature  perturbation  field  plot  at 
t=l .  3  hours.  View  is  zoomed  in  from  left  to 
right,  showing  further  detail  of  the 
convective  dynamic. 

2 .  Flux  Trends  and  Dependencies 

Temperature  flux  averaged  over  the  whole  computational 
domain  was  recorded  and  smoothed  in  time  using  the  running 
average  period  of  20  non-dimensional  units  to  remove  high- 
frequency  oscillations.  Figure  12  shows  the  typical 
evolution  of  heat  flux  as  a  function  of  time.  The 
simulation  is  initiated  at  rest,  followed  by  random  noise 
"kicking"  off  the  double-diffusion  process.  This  is  why 
there  is  strong  flux  growth  almost  immediately  after  the 
simulation  begins,  indicative  of  a  large  amount  of 
potential  energy  being  released  by  the  temperature  field. 
As  time  increases  the  system  reaches  a  state  of  near¬ 
equilibrium,  with  fluxes  fluctuating  about  a  mean  value.  In 
this  thesis,  the  mean  heat  fluxes  were  determined  from  the 
equilibrated  portion  of  the  graph. 
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Heat  Flux  vs.  Time,  (R^=4) 


Figure  12.  A  typical  evolutionary  pattern  of 
temperature  fluxes  as  observed  in 
simulations.  Taken  from  the  experiment  with 
Rp=4. 

Through  one-step  numerical  experiments,  it  is  possible 
to  determine  the  magnitude  of  diffusive  fluxes  and  compare 
them  with  those  derived  from  laboratory  experiments.  Since 
DNS  flux  results  are  non-dimensional,  a  conversion  was  made 
to  get  dimensional  units.  Non-dimensional  and  dimensional 
DNS  temperature  fluxes  are  related  by 

„  dT 

FT  dim  —  ^Tnon-dim^T  ,  (19) 

dz 

with  the  corresponding  (dimensional)  heat  flux  given  by 

FH6im=FT<HmCPP  (20) 

where  Cp  is  the  specific  heat  for  sea  water  and  p  is  the 
density  of  sea  water.  Substituting  (19)  into  (20)  and 
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utilizing  values  calculated  from  ITP  data  for  the  Beaufort 
Gyre  yields  a  conversion  factor  of  0.0075  W/m2  (Prikasky 
2007 ;  Caro  2009) . 

As  indicated  by  Equations  (1)  —  (3) ,  fluxes  are  expected 
to  be  dependent  on  the  density  ratio  for  a  fluid  with  a 
temperature  variation  across  the  interface  (Turner  1973). 
To  determine  this  dependency,  six  one-layer  DNS  were 
performed,  each  with  a  different  Rp  chosen  from 
[2,3,4,6,8,10]  and  HnOndim=300  (~4  m  in  dimensional  units). 
Each  DNS  illustrated  the  same  behavior  as  seen  in  Figure 
12.  As  expected,  fluxes  consistently  decreased  with 
increasing  Rp,  with  measured  values 

F  Hdim=  [0.634,0.527,0.398,0.360,0.316,0.300]  W/rrr. 

To  determine  flux  dependency  on  H,  four  one-layer  DNS 
were  performed  for  Hnondim=  [37 . 5, 75,  150, 300]  ,  equivalent  to 
approximate  dimensional  values  of  Hdim=  [  0 . 5,  1 , 2 , 4  ]  m,  and 
fixed  Rp=4 .  An  increase  in  H  for  fixed  temperature  gradient 
corresponds  to  variation  in  AT,  thus  allowing  systematic 
analysis  of  the  corresponding  pattern  of  variation  of  heat 
transport.  As  anticipated,  fluxes  increased  with  the 
increasing  layer  thickness,  with 

FHdim=  [0 . 053,0.102,0.204,0.398]  W/m2. 

These  dependencies,  and  how  they  compare  to  laboratory 
experiment  results,  are  discussed  next  in  greater  detail, 
in  the  context  of  3D  DNS. 

B.  THREE-DIMENSIONAL  MODEL 

Figure  13  presents  the  typical  distribution  of  the 
total  temperature  (T  +T )  field  realized  in  3D  DNS.  This 
visualization  reflects  the  interaction  between  regular. 
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nearly  flat  interfaces  and  disorganized  convective  plumes 
advecting  temperature  and  salinity  through  well-mixed 
layers.  Shown  are  two  identical  layers  to  illustrate  the 
periodic  nature  of  the  system  (based  on  a  single  layer 
DNS)  . 


Figure  13.  3D  volume  rendering  of  a  diffusive 

convection  temperature  field  for  Rp=4 . 

Graphic  generated  using  VAPOR  software. 

As  was  the  case  with  2D  DNS,  a  series  of  3D 

experiments  were  performed  in  the  configuration  that  is 
consistent  with  typical  staircases  in  the  Beaufort  Gyre  of 
the  Arctic  Ocean  for  various  density  ratios, 

Rp=  [2, 3, 4,  6,  8,  10]  and  step  heights,  Hdim=  [ 0 . 5,  1 , 2 , 4  ]  m.  Refer 
to  Table  1  for  further  details  on  these  nine  experiments. 

As  expected,  heat  flux  values  for  each  DNS  evolved  in  the 
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manner  illustrated  in  Figure  12.  Thus,  the  same  procedure 
for  calculating  mean  fluxes  applied  to  2D  experimentation 
could  be  used  for  all  DNS,  over  similar  time  intervals 
(thereby  ensuring  consistency  in  our  approach) . 

1 .  Flux  Law  Comparisons 

Marmorino  and  Caldwell's  (1976)  and  Kelley's  (1990) 
lab-derived  formulations  of  the  4/3  flux  law  (Equations  (7) 
and  (8),  respectively),  were  used  for  comparison  with  the 
3D  DNS  results.  Figure  14  compares  heat  fluxes  diagnosed 
from  the  DNS  with  those  derived  using  the  4/3  flux  law  for 
a  range  of  density  ratios.  This  comparison  indicates  that 
for  low  Rp,  laboratory  estimates  enjoy  some  success  in 
estimating  heat  flux  (FH)  .  For  instance,  at  Rp=2,  Kelley's 
(1990)  law  is  within  10%  of  the  DNS  value,  with  Marmorino 
and  Caldwell's  (1976)  law  performing  at  a  similar  level  for 
Rp=3 .  However,  for  larger  Rp  (>~3) ,  both  laws  underestimate 
Fh,  with  the  discrepancy  increasing  with  Rp.  An  approximate 
factor  of  2  difference  can  be  seen  at  large  values  of 
density  ratio  (Rp>5) .  This  discrepancy  is  not  surprising, 
given  the  substantial  differences  in  the  numerical  and 
laboratory  setup. 


29 


Heat  Flux  vs.  Density  Ratio 


Figure  14.  A  comparison  of  heat  fluxes  (as  a  function 
of  Rp)  diagnosed  from  3D  DNS  with  those 
derived  from  extant  laboratory-derived  flux 
laws . 

In  order  to  present  the  DNS  findings  in  explicit  form, 
also  presented  in  Figure  14  is  the  best-fit  curve  of  the 
form 

^„o„dhn=«+K^-l)C  (21) 

where  a=31,  b=103.4  and  c=-1.24.  If  we  additionally  assume 

that,  at  the  leading  order,  the  dependence  on  AT  follows  the 
four-thirds  power  law,  then  (21)  leads  to  the  following 
parameterization  of  vertical  transport 

^rd1m=K+^i(^-l)C](^)4/3  (22) 

The  discrepancy  in  the  lab-derived  flux  values  in 


calculating  fluxes  points  to  an  inaccurate  prefactor  value, 
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with  an  error  of  2-3.  Note  that  Equation  (22)  offers  an 
adequate  description  of  the  numerical  data  provided  that 
the  dependence  on  temperature  variation  is  well-described 
by  the  4/3  power  law.  Whether  or  not  this  assumption  is 
accurate  can  also  be  questioned  (Kelley  et  al .  2003) . 

2.  Testing  the  4/3  Exponent 

In  order  to  determine  whether  the  numerical 
simulations  are  adequately  represented  by  the  4/3  power  law 
in  Equation  (3) ,  a  set  of  four  one-step  3D  DNS  with 
different  step-heights  was  performed  with  varying  H  and 
fixed  Rp=4 .  All  other  parameter  settings  were  kept 

identical  for  each.  The  increase  in  H  for  the  fixed 

temperature  gradient  (Tz=l . 361xl0-2  °C/m)  corresponded  to 

the  variation  in  AT.  Figure  15  is  the  plot  of  aFT  as  a 
function  of  aAT  (dimensional  units)  in  logarithmic 
coordinates  as  derived  from  DNS  and  as  estimated  by 
Kelley's  (1990)  law.  From  inspection  of  Figure  15,  it  is 

immediately  apparent  that  the  4/3  exponent  does  not 
accurately  describe  the  numerical  results.  A  best-fit  curve 
of  the  form 

log  (aFT )  =  c  +  m  \og(aAT)  {23) 

yields  m=1.02  (with  c=-13) ,  which  is  significantly  less 
than  the  commonly  assumed  4/3.  If  the  fit  is  made  so  as  to 
exclude  the  lowest  aAT,  then  m=l.ll.  It  also  should  be 
noted  that  laboratory  flux  estimates  came  close  to  the 
numerical  values  for  relatively  large  AT  (equivalent  to 

large  step  heights,  H)  ,  but  differed  for  relatively  thin 
steps . 
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aFTvs.  aAT 


Figure  15.  Testing  "four-thirds  flux  law"  exponent.  aFT 
(m/s)  as  a  function  of  aAT  in  logarithmic 
coordinates,  as  diagnosed  from  3D  DNS. 

C.  COMPARISON  OF  TWO-  AND  THREE-DIMENSIONAL  EXPERIMENTS 

Utilizing  2D  DNS  has  its  obvious  advantages, 
especially  in  a  computationally-restrictive  modeling 
environment.  2D  DNS  are  less  expensive  for  realistic 
parameters,  but  it  is  reasonable  to  question  whether  they 
provide  adequate  approximation  to  their  3D  counterparts. 
For  instance,  it  is  widely-known  (e.g.  Kolmogorov  1941; 
Batchelor  1953)  that  energy  cascades  from  larger  to  smaller 
scales  in  the  full  3D  regime  whereas  the  opposite  is  true 
for  2D  simulations.  To  test  their  accuracy,  2D  DNS  with 
identical  parameter  sets  and  resolutions  to  those  3D  DNS 
previously  described  were  conducted.  Mean  temperature 
fluxes  were  noted  and  compared  by  looking  at  the  resultant 
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3D/2D  ratios  ( FT3D/FT2D )  .  Figure  16  shows  this  ratio  for  the 
DNS  with  varying  Rp  and  fixed  H  (top  panel) ,  and  varying  H 
and  fixed  Rp=4  (bottom  panel)  . 

Surprisingly,  in  both  cases,  2D  DNS  stand  up  quite 
well.  For  varying  Rp,  the  maximum  value  FT3D  /  FT2D  =  1 . 5  7  4 

occurs  at  Rp=2 .  However,  for  all  larger  Rp  values,  2D 
estimates  are  within  10%  of  their  3D  counterparts.  The  same 
consistency  between  2D  and  3D  results  is  seen  when  density 
ratio  is  kept  constant  and  layer  thickness  is  systemically 
varied . 
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Figure  16.  Comparison  of  3D  and  2D  DNS  as  a  function  of 
Rp  (top  panel)  and  H  (bottom  panel) .  With  the 
exception  of  Rp=2,  the  3D/2D  ratios  were 
within  10%  of  FT3D  /  FT2D  =  1  . 
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D.  EFFECTS  OF  BOUNDARY  CONDITIONS 

Carpenter  et  al .  (2012)  performed  a  series  of  3D  DNS 
which  differed  from  those  of  Caro  (2009)  in  that  their 
model  configuration  included  rigid  boundaries.  As  a  result, 
their  DNS  experienced  a  run-down  evolution  of  the  diffusive 
interface,  a  feature  that  is  compatible  with  the  laboratory 
configuration  but  which  may  not  be  realized  in  actual 
oceanic  staircases.  Although  Carpenter  et  al .  (2012)  showed 
that  the  laboratory-based  flux  laws  predict  lower  heat 
transport  than  those  observed  in  their  simulations,  it  is 
likely  that  their  derived  flux  values  are  still  too  low 
because  of  the  implementation  of  rigid  boundary  conditions. 

To  test  the  effects  of  boundaries,  a  numerical 
experiment  with  rigid  stress-free  boundaries  in  all  spatial 
directions  was  conducted;  the  setup  that  reflects  the 
dynamics  of  a  laboratory  run-down  experiment.  The  DNS  were 
run  in  two  dimensions  only,  since  they  have  been  shown  to 
be  generally  adeguate  in  capturing  the  dynamics  of 
diffusive  convection  and  the  associated  heat  transport.  To 
make  the  comparison  between  the  periodic  and  run-down 
simulations  objective,  in  both  cases  the  same  parameters 
were  used  (see  Table  1)  .  At  initialization,  Rp=1.5,  and  as 
the  system  ran  down,  Rp  increased  with  time.  This  made  it 
possible  to  calculate  the  mean  temperature  flux  as  a 
function  of  density  ratio  over  the  continuous  range 
1 . 5<RP<1 0  . 

In  Figure  17  heat  flux  values  as  a  function  of  Rp  from 
this  DNS  are  compared  to  those  derived  from  the  periodic  2D 
simulations.  In  the  run-down  experiments,  heat  flux  was 
computed  at  the  level  of  diffusive  interface  (z=H/2)  and 
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smoothed  in  time  using  the  running  average  interval  of  100. 
For  smaller  values  of  Rp  (<3) ,  fluxes  with  rigid  boundaries 
are  larger  with  the  difference  rapidly  diminishing  as  Rp 
increases.  Quickly,  the  run-down  fluxes  fall  significantly 
below  the  corresponding  estimates  from  the  periodic 
simulations.  For  comparison  purposes,  flux  values  as 
determined  by  Kelley's  (1990)  law  were  included.  As  is  to 
be  expected  from  a  flux-law  obtained  in  the  laboratory 
setting  (where  rigid  boundaries  are  present),  Kelley's  law 
provides  better  overall  approximation  to  the  run-down  DNS- 
based  flux  values,  than  the  periodic  system. 


Heat  Flux  vs.  Density  Ratio 


Figure  17.  The  effect  of  boundary  conditions  as  seen 
when  comparing  fluxes  as  a  function  of  Rp.  As 
Rp  increases,  fluxes  without  boundaries  are 
significantly  higher  than  their  non-periodic 
counterparts . 
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This  DNS  demonstrates  the  inherent  difficulties  that 
arise  when  studying  diffusive  convection  in  a  laboratory 
setting.  The  rigid  boundaries  and/or  run-down  character  of 
experiments  act  to  significantly  suppress  the  vertical  heat 
transport.  Since  a  rigid  bottom  and  side-walls  are  present 
in  laboratory  experiments,  the  application  of  the 

laboratory-derived  flux  laws  to  oceanic  interfaces,  which 
are  not  constrained  by  rigid  boundaries  and  which  can 

persist  for  years,  may  not  be  appropriate.  Based  on  the 
high-resolution  experiments  with  periodic  boundary 

conditions  it  is  estimated  that  the  error  of  extant  flux 
laws  is  modest  (a  factor  of  two  or  less)  and  could  be 


tolerated 

in 

many 

applications . 

However, 

if  a 

more  accurate 

estimate 

of 

the 

vertical  heat 

transport 

is 

needed,  using 

Equation 

(22 

)  is 

recommended. 

since  it 

is 

based  on  the 

periodic  simulation  results. 
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IV.  DATA  AND  DATA  PROCESSING 


A.  ICE-TETHERED  PROFILERS  AND  PREVIOUS  RESEARCH 

The  Woods  Hole  Oceanographic  Institute  (WHOI)  ITP 
program  provides  timely  and  accurate  oceanographic 
measurements  from  the  Arctic  Ocean.  ITPs  assist  in  assuring 
consistent  measuring  capability  in  regions  that  have,  until 
now,  remained  data-sparse  due  to  access  limitations.  The 
ITP  is  a  moored  conductivity,  temperature,  and  depth  (CTD) 
sensor  that  returns  high  vertical  resolution  measurements 
of  upper  ocean  temperature  and  salinity  beneath  sea-ice. 
Each  one  is  designed  to  move  with  the  ice  it  is  moored  to 
and  consists  of  a  small  platform  on  top  of  which  rests  a 
buoy.  A  10-inch  hole  is  drilled  through  the  sea  ice  and  the 
instrument  is  inserted  in  the  water  from  above.  The 
instrumentation  is  tethered  to  the  surface  float  via  a  wire 
rope  that  is  weighted  on  the  bottom.  This  rope  is 
approximately  800  meters  in  length,  and  the  instrument 
traverses  the  wire  rope  several  times  daily  (Figure  18)  . 
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Figure  18.  Ice-Tethered  Profiler  (ITP)  schematic  (WHOI 
ITP  website  2012)  . 


Data  from  ITPs  1-6  were  used  for  comparison  with 
previous  studies  that  focused  on  the  double-diffusive 
staircases  in  the  Beaufort  Gyre  (e.g.  Timmermans  et  al . 
2008;  Caro  2009;  Chaplin  2009)  (see  Figure  19) .  Of  specific 
interest,  Timmermans  et  al .  (2008)  used  the  potential 
temperature  and  salinity  measurements  recorded  by  the  six 
ITPs  to  estimate  vertical  heat  fluxes  via  Kelley's  (1990) 
laboratory-derived  flux  law.  Figures  20  and  21  present  the 
maps  of  average  density  ratio  and  associated  heat  flux  as 
calculated  in  Timmermans  et  al .  (2008)  . 

Their  study  concluded  that  the  vertical  transport  of 

heat  from  the  AW  is  unlikely  to  have  a  significant  impact 
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on  the  Arctic  Ocean  surface  heat  budget.  Estimated  vertical 
heat  fluxes  were  found  to  be  in  the  range  0.05-0.3  W/m2, 
with  the  average  value  of  0.22  W/m2;  this  is  considerably 
less  than  the  values  derived  from  our  3D  DNS  results. 


Figure  19.  Trajectories  of  ITPs  1-6  (WHOI  ITP  website 
2012)  . 


By  analyzing  the  same  data  as  Timmermans  et  al .  (2008) 
and  applying  the  "new"  flux  law  (22),  a  similar,  but  more 
up-to-date  counterpart  to  Figure  21  can  be  obtained. 
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Figure  20.  Rp  averaged  over  the  thermocline  contoured 
from  ITP  measurements  (After  Timmermans  et 
al .  2008)  . 


Figure  21.  Map  of  heat  flux  (W/m2)  estimated  by  ITPs  1-6 
profiles,  averaged  over  the  200  m  to  300  m 
deep  thermohaline  staircase  (Timmermans  et 
al.  2008)  . 
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B.  DATA  PROCESSING 

The  ITP  data  was  arranged  by  profiler  number  1-6  and 
processed  together  as  one  data  set.  Table  2  outlines  the 
properties  of  these  six  ITPs.  Processing  began  with  the 
corrected  ITP  data  obtained  from  the  WHOI  ITP  website. 
Wilson  (2007)  suggested  that  anomalous  results  were 
obtained  from  the  profiles  corresponding  to  the  downward 
data  collection  for  each  cast.  Possibly  due  to  sensor 
contamination,  these  measurements  were  affected,  and  all 
the  downward  collected  profiles  were  removed  from  the  total 
dataset.  Figure  22  displays  the  location  for  each  profile 
as  well  as  the  data-dense  region  chosen  for  analysis. 


ITP 

#  of  profiles 

Duration 

Distance 

Launch  Date 

Recovery  Date 

1 

2043 

723 

4990 

8/15/05 

8/8/07 

2 

244 

40 

378 

8/19/04 

9/28/04 

3 

1532 

382 

2559 

8/23/05 

9/9/06 

4 

698 

349 

2485 

9/3/06 

8/17/07 

5 

1095 

450 

3892 

9/7/06 

12/1/07 

6 

1335 

1163 

7353 

9/4/06 

10/11/09 

Table  2.  Properties  of  ITPs  1-6. 
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Figure  22.  Locations  for  all  ITP  profiles.  Red 

asterisks  indicate  the  area  analyzed  in  this 
study  using  the  Matlab  scattered  data 
interpolation  function. 

The  majority  of  diffusive  layers  occur  in  the  200-300 
m  region  of  the  thermocline.  The  dataset  was  then  reduced 
to  only  those  data  from  200-300  db  pressures,  which 
corresponds  to  roughly  195-295  m  depth.  By  finding  the 
local  maximum  value  of  potential  temperature  for  each 
layer,  the  center  of  each  layer  was  determined.  Next  the 
diffusive  layers  were  classified  independently  of  the 
particular  cast  using  T-S  plots.  These  plots  (Figure  23) 
show  vertical  clusters  at  nearly  constant  salinity  values 
that  are  indicative  of  horizontally  homogeneous,  diffusive 
layers.  Due  to  the  remarkable  lateral  coherence  of  layers 
in  the  thermohaline  staircase,  it  was  possible  to 
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consolidate  all  six  profilers'  output  into  one  dataset,  by 
identifying  the  pronounced  layers  that  appear  uniformly 
throughout  all  the  data  (Chaplin  2009) . 


Figure  23.  Temperature— Salinity  plots  for  ITPs  1-6 
(After  Chaplin  2009) . 

This  consolidation  of  data  revealed  a  total  of  19 
vertical  layers  common  to  each  profiler.  Using  the 
temperature  and  salinity  jumps  across  each  step,  heat 
fluxes  and  density  ratios  were  computed  for  each  layer, 
and  then  vertically  averaged  over  the  entire  thermocline. 
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c. 


RESULTS 


Heat  fluxes  were  derived  using  Kelley's  (1990)  flux 
law  (8),  and  compared  to  those  computed  using  Equation 

(22) .  Figure  24  is  a  contour  map  of  R  that  can  be  used  to 

quickly  draw  parallels  between  areas  of  low  (high)  density 
ratio  and  high  (low)  heat  fluxes.  Figures  25  and  26  are 
contour  maps  of  heat  flux  colored  using  an  identical  scale 
for  simple  comparison  between  the  two  flux  laws.  Kelley's 
(1990)  law  used  in  Figure  25  under-predicts  heat  fluxes,  on 
average,  by  a  factor  of  two. 

Density  Ratio  {Rp) 
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6 
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Longitude 

Figure  24.  Contour  map  of  R  averaged  over  the  19  layers 
detected  in  the  thermohaline  staircase. 
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Figure  25.  Contour  map  of  heat  flux  (W/m2)  as  computed 
by  Kelley's  (1990)  four-thirds  flux  law, 
averaged  over  the  19  layers  detected  in  the 
thermohaline  staircase 
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Figure  26.  Contour  map  of  heat  flux  (W/m2)  as  computed 
by  Equation  (22),  averaged  over  the  19 
layers  detected  in  the  thermohaline 
staircase 
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V.  CONCLUSIONS 


A.  DISCUSSION 

3D  DNS  are  proving  to  be  the  new  standard  for  analysis 
of  thermohaline  staircases  and  the  advective  processes  that 
transport  heat  upwards  across  diffusive  interfaces. 
However,  in  order  to  provide  reliable  estimates  of  heat 

transport,  simulations  should  meet  a  set  of  stringent 
requirements.  In  particular,  simulations  should 

•  Be  performed  at  a  resolution  that  adequately 

represents  the  scales  of  salt  dissipation 

•  Utilize  sufficiently  large  computational  domains 

which  replicate  typical  step-sizes  in  the  Arctic 
staircase 

•  Be  performed  over  a  time  period  long  enough  to 

provide  accurate  statistical  description 

•  Use  parameters  and  conditions,  such  as  a  lack  of 
rigid  boundaries,  reflective  of  those  found  in 
the  Arctic  thermocline 

Such  requirements  are  so  computationally-demanding 
that,  until  now,  DNS  of  Arctic  staircases  have  been 
unfeasible.  The  fulfillment  of  such  rigorous  standards  has 
allowed  several  important  conclusions  to  be  drawn. 

It  was  determined  that  the  extant  laboratory-derived 
flux  laws,  while  qualitatively  consistent  with  numerical 
inferences,  underestimate  the  vertical  heat  transport  by  as 
much  as  a  factor  of  two.  An  alternative  DNS-based  flux  law 
was  formulated.  The  exponent  of  Turner's  (1965)  four-thirds 
flux  law  was  also  tested  and  found  to  not  offer  an  optimal 
description  of  the  numerical  data.  Evidence  suggests  that 
the  4/3  exponent  in  Equation  (3)  may  require  downward 
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adjustment  by  as  much  as  30%,  if  it  is  to  be  used  for  flux 
calculations . 

Realistic  3D  simulations  which  are  inherently 
computationally  expensive,  were  found  to  be  consistent  with 
the  equivalent  2D  DNS.  The  latter  can  provide  an  attractive 
alternative  for  sufficiently  large  density  ratios  (Rph3) 
but  the  margin  of  error  increases  with  smaller  density 
ratio  values  (Rp<2) . 

The  choice  of  boundary  conditions  (rigid  versus 
periodic)  and  the  experimental  setup  (run-down  versus 
maintained)  had  a  substantial  effect  on  the  vertical  heat 
transport.  It  is  suggested  that  oceanic  staircases  are 
better  represented  by  models  with  maintained  periodic 
boundary  conditions. 

Finally,  Equation  (22)  was  applied  to  ITP  data  in  the 
vicinity  of  the  Beaufort  Gyre  and  the  results  compared  with 
those  from  previous  heat  flux  studies.  It  has  been  shown 
that  the  four-thirds  flux  law  calculations  underestimate 
heat  flux  by  a  factor  of  two  across  a  range  of  density 
ratios  representative  of  the  Arctic  thermocline,  with 
Fh=0.4  4  W/m2. 

This  study  was  prompted  by  the  melting  of  Arctic  sea- 
ice  and  the  question  as  to  whether  double-diffusive 
convection  can  adequately  shuttle  AW  water  vertically  to 
impact  the  melt  rate  of  surface  ice.  Although  the  entire 
surplus  of  1  W/m2  of  heat  in  the  Arctic  heat  budget  is  not 
fully  accounted  for,  our  findings  do  help  to  close  the  gap. 
That  vertical  heat  flux  is  (on  average)  twice  as  large  as 
previously  believed,  points  to  diffusive  convection  playing 

a  larger  role  in  transporting  AW  vertically  to  the  surface. 
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Alone  or  in  combination,  taking  into  account  all  ocean  and 
atmospheric  processes  could  explain  the  existing  ice-melt 
observations  and  heat  budget  surplus.  Furthermore,  if  AW 
continues  to  warm,  double-diffusive  convection  will  likely 
play  a  larger  role  in  increasing  the  rate  at  which  surface 
ice  melts. 

B .  RECOMMENDATIONS 

The  analysis  in  this  thesis  warrants  future  work  on 
the  subject.  Through  further  application  of  the  DNS-derived 
flux  law  (22)  to  ITP  data,  the  scope  of  study  can  be 
increased.  For  instance,  the  analysis  of  more  recent  data 
profiles  can  help  determine  the  role  of  a  continually- 
warming  AW  on  the  Arctic  sea-ice  melt  rate.  By  taking  a 
mean  average,  it  would  be  relatively  simple  to  determine 
the  changes  in  vertical  heat  transport  through  Arctic 
thermohaline  staircases  over  time. 

Re-engaging  in  the  inverse  modeling  technique  that 
was  started  by  Chaplin  (2009)  would  also  help  to  infer  the 
vertical  heat  flux  values  through  diffusive  staircases. 
Previous  findings  support  the  DNS  results  from  this  study, 
but  also  raised  questions  on  the  algorithm  performance  and 
the  sensitivities  associated  with  the  models  due  to  varying 
layer  thickness.  These  issues  should  be  addressed  in  future 
research . 

Finally,  the  newly  derived  flux  parameterizations 
should  be  considered  for  inclusion  into  lower-resolution 
basin  scale  simulations  in  order  to  observe  the  impact  of 
double  diffusion  on  large-scale  circulation  patterns. 
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